clc;
clear;
format short

%% 网格离散数据
JM = [1 2 3; 3 4 5; 5 6 7];
JX = (0:1/6:1)';
JB1 = [1, 1];
JB2 = [3, 3, 2];
N = 7; % 节点总数
E = 3; % 单元总数

%% K 和 b 的计算
K = zeros(N, N); % 初始化 K
b = zeros(N, 1); % 初始化 b
for i = 1:E
    for j = 1:3 % 取出第 i 个单元三个节点的坐标
        jiedian_haoma(j) = JM(i, j);
        JXe(j, 1) = JX(jiedian_haoma(j), 1);
    end
    ke = function_of_ke(JXe); % 调用 function_of_ke 函数, 计算单元 ke
    be = function_of_be(JXe); % 调用 function_of_be 函数, 计算单元 be
    for m = 1:3 % 双层循环组合 K
        for n = 1:3
            K(JM(i, m), JM(i, n)) = K(JM(i, m), JM(i, n)) + ke(m, n);
        end
    end
    for m = 1:3 % 单层循环组合 b
        b(JM(i, m), 1) = b(JM(i, m), 1) + be(m, 1);
    end
end

%% 代入 JB2, 计算 F
F = zeros(7, 1); % 初始化 F
for i = 1:length(JB2(:, 1)) % 循环次数为 JB2 的行数
    II = JB2(i, 1); % 提取边界单元编号
    Fe = function_of_Fe(JB2(i, :)); % 调用 function_of_Fe 函数, 计算单元 Fe
    for m = 1:3 % 根据 JM, 将 Fe 组合到 F 中
        F(JM(II, m), 1) = F(JM(II, m), 1) + Fe(m, 1); % 式(1.74)
    end
end

%% b 和 F 相加
bF = b + F;

%% 采用乘大数带入法, 代入 JB1
for i = 1:length(JB1(:, 1)) % 循环次数为 JB1 的行数
    II = JB1(i, 1); % 提取边界节点编号
    K(II, II) = K(II, II)*1e16; % 乘大数代入法第一步
    bF(II, 1) = K(II, II)*JB1(i, 2); % 乘大数代入法第一步
end

%% 求解方程, 结果对比
u = K\bF; % 方程求解
result = [JX, u] % 结果输出
plot(JX, u, '*'); % 数值解绘图
hold on % 保持绘图框为打开状态
x = 0:0.01:1; % 精确解定义域
y = 1/2*x.^2 + x + 1; % 精确解表达式
plot(x, y, 'r'); % 精确解绘图

